Local versus global equilibration near the bosonic Mott-superfluid transition 



Stefan S. Natu, 1 ^ Kaden R. A. Hazzard, 2 and Erich J. Mueller 1 

Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853, USA. 
2 JILA and University of Colorado, Boulder, Colorado 80309-0440, USA 

We study the response of trapped two dimensional cold bosons to time dependent lattices. We 
find that in lattice ramps from 11 (superfluid, h/Ui = 3ms, h/Ji = 45ms) to 16 recoils (Mott, 
h/Uf = 2ms, h/Jf = 130ms) the local number fluctuations remains at their equilibrium values if 
ramps are slower than 3 ms. Global transport, however, is much slower (Is), especially in the 
presence of Mott shells. This separation of timescales has practical implications for cold atom 
experiments and cooling protocols. 



Introduction. — Understanding and controlling the 
equilibration of cold atom systems is one of the most 
important current challenges in the field. As isolated 
systems, the relaxation mechanisms are intrinsic and fun- 
damental 0-0]. The atomic systems are readily driven 
far from equilibrium, and are well suited for quantifying 
concepts of non-equilibrium dynamics Moreover, 
controlling the equilibration of cold atoms is key to the 
next generation of experiments: for example one needs 
fast equilibration for condensed matter emulators [l3j |. 
Motivated by recent experiments 14H16I] , we conduct nu- 
merical simulations of the response of a gas of bosons to 
a change in the intensity of an applied optical lattice. 

Despite being performed under similar conditions, 
three recent experiments 14H16I] find relaxation rates for 
two-dimensional lattice bosons that differ by two orders 
of magnitude. Here we show that these discrepancies 
can be explained by a separation of timescales for local 
equilibration and global transport. We illustrate this re- 
sult by numerical simulations within a time-dependent 
Gutzwiller mean-field theory. We further explore the pa- 
rameters, such as system size and trap geometry, which 
influence these timescales. 

The separation of timescales for local and global equi- 
librium is unsurprising, and emerges in most interacting 
systems and materials. For example, in the air around us, 
local equilibrium is achieved on the collision time (~ns) , 
but global equilibrium is limited by transport coefficients 
and is relatively slow. Typically one expects the slow 
variables to be those that are conserved (such as den- 
sity and energy density) and those which correspond to 
broken symmetries (such as the phase of the superfluid 
order parameter). Although we do not do so here, inte- 
grating out the fast degrees of freedom leaves "hydrody- 
namic" equations for the slow degrees of freedom. The 
form of these hydrodynamic equations are strongly con- 
strained b y symm etries, allowing phenomenological de- 
scriptions 



17, 1 



A practical consequence of this separation of timescales 
is that adiabaticity is much easier to maintain if one 
changes parameters in such a way that very little mass 
transport is necessary - a principle which is widely used 
in cold atom experiments. 



Theoretical Setup. — Bosonic atoms trapped by inter- 
fering laser beams are well described by the Bosc Hub- 
bard Hamiltonian UM 
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where a and are bosonic annihilation and creation 
operators, J is the tunneling, and U is on-site interac- 
tion. We denote [i l = /i — V cx (i), where /i is the chemi- 
cal potential and V ex (i) is the external potential at site 
i [20l | . The first sum is over all nearest neighbor sites 
in the plane. In Figure [TJ we show U and J as a func- 
tion of lattice depth V R for 87 Rb in a d = 680nm lat- 
tice generated by light of wavelength A = 1360nm. For 
deep lattices, U = ^/8/^(ka s )E R ^/V R /E R {V Rz /E R ) 1 ^ 

and J = (4/^)(S fl Vi) 1 / 4 exp(-2/lV^), where 

E R = k 2 /2m is the lattice recoil energy in terms of the 
light wave- vector k = 2tt/X for light with wavelength A, 
Vr, V Rz are the radial and axial lattice depths, a s is the 
scattering length j2l|. Different two-dimensional experi- 
ments use different strengths of axial confinement (V Rz ). 
Since U only depends on V^J Z 4 , we will make the simplest 
choice, V Rz = V R . None of our conclusions are qualita- 
tively affected by this assumption. 

We calculate dynamics using a time dependent 
Gutzwiller ansatz |19j . which approximates the wave- 
function by ^ = 0, t J2 m c ™ ( f )l TO )i where \m) l is the 
m-particle Fock state on site i, and the coefficients (t) 
are generally space and time dependent. 

This mean-field ansatz reduces Eq.[T]to a sum of single 
site Hamiltonians Hi = — 4t((a i )*a i + (a l )al)+4:t\{ai)\ 2 + 
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V(i)]rii at each site i. Truncating the 



basis at each site to a maximum M particles, Hi is an 
(M + 1) x (M + 1) matrix at each site, and depends on 
the other sites only through (a*) = (1/4) ^ 



(j)\ a i): where 



(cij) = J2 m Vm+ lc„ +1 Cm , and the sum over j includes 
all four nearest neighbor sites. 

Schrodinger's equation idtip = Hip for $ yields a set 
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FIG. 1: Energy scales as a function of lattice depth: 

Microscopic parameters in the 2D Bose-Hubbard Hamiltonian 
(Eq[T]): 4 J (solid), and U (dashed) as a function of lattice 
depth [3 for 87 Rb in a d = 680 nm lattice. The dotted curves 
are the two lowest k — excitations from linearizing Eq. Q 
at unity filling. In the superfiuid state, the Goldstone mode 
has zero energy. In the Mott state, these are the particle/hole 
excitations. 



of differential equations for the d m : 



*dtc l m {t) = ~4J(t)((a t )Wm + lc l m+1 + {o^y/m^j) + 



m(m-l)-» i m + U(t)\(<*i)\ 2 )'LQ) 



The tunnelings and on-site interactions are dynamically 
tuned by changing the lattice depth in time (t). We study 
population dynamics across the superfluid-insulator tran- 
sition by ramping the lattice linearly in time using the 
protocol V(t) =Vi + (V { - Vi)(t/T T ), where V{ and V f are 
the initial and final lattice depths, and r r is the ramp 
time. We consider a time independent radially symmet- 
ric harmonic trap, V cx = \rrvJ 1 {x 2 + y 2 ). 

We approximate the ground state by finding the sta- 
tionary solution to Eq. @, c l m (t) = e~ l<Lt c l m , where e can 
be identified with the energy per site. We use an itera- 
tive algorithm, starting with a trial on, then find c l m by 
solving the eigenvalue problem in Eq. ([2]). We calculate 
a new a% and repeat until the subsequent change in on 
is sufficiently small. To calculate time dynamics, we use 
a split-step method [22| and sequential site updates 23 1 . 
This approach conserves both total particle number and 
energy (for time- independent Hamiltonians) . 

The resulting dynamics describe the behavior of a sin- 
gle quantum state, rather than a density matrix. Nev- 
ertheless, the equations governing the time dependent 
Gutzwiller ansatz are highly nonlinear and contain a 
large number of degrees of freedom. This structure is 
rich enough that under appropriate conditions time dy- 
namics leads to thermalization, with (on average) energy 
equally distributed among all modes. 

Results. — We consider several different scenarios in or- 
der to fully explore the response this system to a lattice 
ramp. We start by analyzing a homogeneous system: 
this investigation yields the timescale for local equilibra- 
tion. This timescale sets the fundamental limit for how 
fast equilibration can take place in the absence of global 



mass transport. Similar to the Harvard experiments |15| . 
we find that local equilibrium can be maintained even 
under relatively rapid quenches through the superfluid- 
Mott boundary. 

Next we explore the requirements for maintaining 
global equilibrium. We show that equilibration times 
are much longer in systems requiring large amounts of 
particle transport. This situation is exacerbated by the 
presence of large Mott domains, as in the Chicago exper- 
iments 0. 

We conclude by showing that even in large systems, 
rapid global equilibration can be achieved, if the trap pa- 
rameters are chosen in a way as to minimize transport 
between intervening Mott shells. Our results in this sec- 
tion are consistent with the Munich experiments fl6l ] . 

Local equilibration. — In an isolated homogeneous sys- 
tem, ramping the depth of an optical lattice does not 
lead to bulk mass transport. Instead, all of the temporal 
dynamics simply involve the evolution of number fluctua- 
tions and correlations. Thus equilibration is governed by 
local physics and Eq.[2]reduces to the single site problem. 
We numerically integrate this nonlinear set of ordinary 
differential equations, taking J and U functions of time, 
corresponding to a linear ramp of the lattice from depth 
V[ to Vf. We vary V\, Vf, and the ramp time r r . We take 
all parameters to correspond to 87 Rb atoms, and take 
n = 1 particles per site. 

At unity filling, near the Mott transition, we truncate 
the basis to at most 2 particles per site. In this trun- 
cated basis, the probability of having a single particle 
per site P{1) is identical to the probability of having an 
odd number of particles per site, which is the experimen- 
tal observable in the Harvard experiments [lEj . 

Both of the gapped q = single-particle excitations 
(see FigO] and Ref. 24]), and the continuum of two- 
phonon excitations contribute to the non-adiabatic evolu- 
tion. All of these modes are captured in a time-dependent 
Gutzwiller framework 
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One expects that the num- 
ber of excitations goes to zero as the ramp rate van- 
ishes. When gapped excitations of energy A dominate 
the dissipation, then the ramp becomes adiabatic when 
■^dA/dt < 1 (H. 

In Fig. [2] we show that the timescale for local equi- 
libration is very short. Starting with a superfiuid at 
Vi = 11-Er, we ramp up to different lattice depths. We 
plot the time evolution of the probability that a single 
particle sits at a given site as we vary the the ramp time r r 
from 0.1h/Ui to lOfi/t/j, where Z7j = h/3ms. Our scheme 
is identical to that considered by the Harvard experi- 
ments [l5|. Fitting these curves to simple exponentials 
yields a characteristic timescale r a , which as we show in 
the inset, is comparable to . 

Inhomogeneous dynamics. — We now consider an inho- 
mogenous system by imposing a harmonic external po- 
tential on top of the lattice. The protocol for lattice 
ramps is same as before, starting with a superfiuid at 
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FIG. 2: (Color Online) Population dynamics at unity 
density n = 1 (Top): Probability of having one particle per 
site at the end of a lattice ramp from Vi = UEr lattice to (top 
to bottom) V{ — 13(yellow), 15(green), 17(blue), 19(purple) 
and 25 (red) in units of Er after different lattice ramp times 
r r = 0.1/Ui ~ 0.3ms to 10/f/i. Inset: Fitting these curves to 
simple exponentials yields a fast timescale for lattice equili- 
bration of r a ~ The best fit line is shown as a guide to 
the eye. 



11 recoil lattice depth. The central chemical potential is 
chosen such that the central density is close to unity, jus- 
tifying the truncated basis (M = 2) used here. Through- 
out we define time in units of 2ir/Ui where U\ is the on-site 
interaction at V\ = 11E R equal to ~ 2ir x 300Hz. We use 
a trapping frequency u = 15Hz. 

In Fig. |3] we plot the density profile after a lattice ramp 
from Vi = 11E R to V t = 16E R in a time t = 120 x 
2n/Ui for a system 30 x 30 sites containing 500 particles. 
As shown previously, this ramp is sufficiently slow to be 
locally adiabatic. The parameters are chosen such that 
at later times, a large Mott region separates the central 
superfluid from the edge in the final state. In the Chicago 
experiment, [3], this Mott domain was ~ 50 sites wide 
while in our simulations it is ~ 8 sites. 

Like the experiment, we find that after this time the 
density profile of the final state (dashed line) is very dif- 
ferent from the equilibrium state at Vf (dotted red), im- 
plying a relaxation time much longer the ramp time of 
400ms. Indeed, further simulations show that it is longer 
than the experimental timescale of seconds. In the re- 
mainder of this section we describe the cause of the slow 
equilibration, and conduct a number of additional sim- 
ulations to illustrate how equilbration times depend on 
the various experimental parameters. 

The major bottleneck for equilibration in Fig. [3]is mass 
transport across the Mott region. To illustrate the spatial 
location of the Mott insulator, in Fig.|3jb) we plot the co- 
herences Ci = —(di) J2j( a j) as a function of time, where 
i,j denote nearest neighbor pairs. Mott regions (C = 0) 
show up as dark regions in the density plot. The Mott 
plateau widens over time, isolating the central superfluid. 
The peak atomic density in the initial lattice exceeds that 
of the equilibrium state at the final lattice depth. How- 
ever the Mott region prevents mass flow from the center 
to the edge. The exponential suppression of transport 
through the Mott region was studied by Vishveshwara 
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FIG. 3: (Color Online) Slow transport across Mott re- 
gion (Left) Evolution of an initial superfluid state (solid) at 
Vi = HE a and N = 500 in a 15Hz radial trapping potential. 
Final density profile (dashed) afterarampr r = 120x27r/?7i ~ 
400 ms, is very different from the equilibrium state (dotted) 
at V{ — 16Er. (Right) Density plot showing the time evolu- 



tion of the coherences (d 



i a i) J2j{ a *j))> a growing Mott 



region in the wings which cuts off transport in the interven- 
ing superfluid producing a non-equilibrium final state at late 
times. Brighter colors correspond to higher coherence 



and Lannert [26f . 

Fast equilibration without transport. — Here we show 
that equilibration times can be dramatically reduced 
when parameters are chosen such that no bulk transport 
across Mott regions is required. The parameters are cho- 
sen to mimic the large systems considered by Sherson et 
al. (lfj| . which attained global equilibrium on timescales 
comparable to 100ms. Figure [4] shows the time-evolution 
of an initial state at V = 11 at N = 800 in 2.5Hz trap, 
and a central chemical potential of 1AU. We find that 
after an evolution of r r = 25 x 2ir/U- u the final profile 
(dashed) is close to the equilibrium T = Gutzwillcr 
prediction (dotted). 

Despite the fact that the n — 1 Mott region is of similar 
size, we find faster equilibration times in this system as 
compared to the one in Fig|3l The difference is that 
here parameters are chosen such that the total number 
of particles in the center is the same in the initial and 
final states. Thus no transport is needed across the Mott 
reion. 

Summary. — Our work was motivated by three exper- 
iments. The Chicago experiments 14j have a large sys- 
tem, and an extremely wide Mott region (~ 50 sites) 
inhibiting transport. We showed by smaller scale calcu- 
lations that once the Mott shell is a significant fraction 
of the system size, the dynamics slow dramatically, lead- 
ing to extremely long equilibration times observed in the 
experiment. 

The Harvard experiments [15j have good evidence of 
local equilibrium on short timescales. We investigate this 
timescale by studying a homogenous system. One insight 
into the timescale for equilibration is that the lowest en- 
ergy k — single particle excitation has a gap ~ U. This 
energy scale appears to set the timescale for adiabaticity. 
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FIG. 4: (Color Online) Time-evolution at higher den- 

sity(Left): Evolution of an initial superfiuid state for V\ = 
ll-Eij and N = 900 (solid) in a 15Hz radial trapping poten- 
tial in a linear ramp with r r = 25 x 2n/Ui = 80ms. The 
dotted profile is the T = equilibrium Gutzwiller profile at 
Vo = 16Er for the same parameters. The final density profile 
(dashed) agrees with the T = equilibrium Gutzwiller profile. 
(Right) Time evolution of the spatial coherence distribution, 
showing the formation of an n = 1 and n — 2 Mott plateaus. 
Lighter colors imply larger coherences. 

The Munich experiments [16| find excellent agreement 
with equilibrium profiles after very short 75ms (~ J) 
ramps. For the parameters we consider, our simula- 
tions reproduce this result. We attribute the differ- 
ence between the Chicago and Munich observations to 
the greater amount of transport accross Mott regions 
required to reach equilibrium for the Chicago parame- 
ters. We believe that although the Harvard experiments 
have good evidence for local equilibration on timescales 
of 1/U, they may be out of global equilibrium after such 
short times. 
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